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ABSTRACT 

Here is a simple example of control of inherently unstable 
system. An inverted pendulum pivoted on top of a cart is to be 
stabilized by applying force to the cart through an electric motor. 

In the electrical laboratory of the United States Naval Post- 
graduate School, a cart with a stick pivoted on top of it has been 
built, tested and simulated with CDC 1604 digital computer. 

The author, Lieutenant Mu-yu Wan of the Chinese Navy, wishes 
to thank Dr. Harold A. Titus of the United States Naval Postgrad- 
uate School for his patient assistance in this work as thesis super- 


visor. 
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1. Introduction 

It seems to be interesting to provide artificial stability for an 
inherently unstable physical system. Some immediate questions 
arise, suchas: what are the largest errors which can possibly be 
correctéd with a limited available control force, and what is the 
best control strategy which minimizes the time required and maxi- 
mizes the size of the system errors which can be corrected. 

In this simple case here, an inverted pendulum pivoted on top 
of a cart is to be stabilized by applying limited force to the cart 
through an electric motor. This inherently unstable system is 
assumed to be adequately represented by a set of linear differen- 
tial equations over the range of interest. The type of instability 
is represented by real non-negative characteristic roots. The 
motor supply voltage is manipulated, within its allowable limits, 
as a function of the state variables of the set of linear differential 
equations. 

In Chapter 3, an analog computer is used to really balance 
the pendulum. The general schematic of the system is shown here. 


control [unstable — = 
Computer system m 


state 
፡ variables 










force 






- 


FIG. (1-1). General Schematic 
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In Chapter 4, the whole system is simulated with CDC 1604 
digital computer, Graphs are plotted and situations discussed. 

In Chapter 5, the technique;of optimal discrete-time control 
is introduced. The results simulated by CDC 1604 turned out to 


be successful. 





2. Uncontrolled System 
2. 1 Linear Differential Equations 

As shown in Fig (2-1-1), the instantaneous angle that the pendu- 
lum makes with its unstable equilibrium position is 0, and the 
position of the cart with respect to some reference point on the 


floor is Е... The coordinates are shown with positive displacement. 


-----5 + moto yoltage V 


FIG. (2-1-1). System To Be Stabilized 
For establishing the equations of motion, we define some addi- 
tional symbols: 
f force applied to cart 
ውን damping coefficient including friction and back e. m. f. 


fv applied voltage force coefficient 





g acceleration of gravity 
M total system effective mass 
m mass of the pendulum 
r distance of pendulum mass center from hinge line 
y pendulum radius of gyration about hinge line 
V  vpltage applied to d. c. drive motor 
The equations of motion can be found by consideration of the 


following diagram. 





FIG. (2-1-2). Force Diagram Showing Equations of Motion 
* Note: 8 is very small. 


The liniarized equations of motion are 


2. 4. 
mı O0 mIEd rs (271-1) 





and 


M É = - mrŠ + f CDL) 
where 
ст с Ga) 


We have data measurements as follows: 


m 
pa (а ee 


npe ovs 

N = 6:10 kg. ' 

r = 0,44 m. 

l = length of longer part from hinge = 0. 94 m. 


Be length oi shorter part from hinge UO mi. 
2. 

A A 0 254 

V = 24 volts (for selected D. C. motor) 


The force exerted on the cart is 13. 6 newton while the cart is 


held still, i.e., 5 = 0, thus from (2. 1. 3) 


In (2.1.2), we have mass of the pendulum much less than mass of 
the whole system, the term "mr" can be neglected, given: 
ыд (2. 1.4) 
Securing the pendulum (8 = 0), we run the cart on the floor 
and observe the motion with a brush recorder, find the average 


velocity and acceleration as: 





© = 0.90 ™/sec 


and 
с zn m/ sec 
By (2. 1. 4) 
a E Ес повест 
š : 
Now (2. I. 1) and (2.1. 4) can be written as: 
, А Y4 Е Su T. 
9 y? ፀ 2 Е ; 4 (2.1.1) 
и 163% “መ ጋ 
Е = 7 E V ( ) 
Define 
0-80 1 
0 = 91 = Ө» 
=, 
EE 
Then, we get à set of linear differential equations as: 
9 = 35 
TOE 
& m6 Rpg, a V 
Ў = 52 
2 => 1155)" AV] 
In matrix form: 
9) መ ው) 9, 0 
XM Th EU 
99 $2 0 0 Me: 9, E $: E OS 
ООО ғ, o | Mg 





Substituting numberical values: 


91 . ЕЕЕ: 9, 0 

ፀ 72 0 0 3.5 ፀ SUD (2.1.5) 
E 5927) E, ЕП 

4 | pL 0-2 E, 9. 8 


4 
Or, by defining some new matrix names, and х! s as the name of 


state variables: 


x= [A] x te u (2.1. 6) 
where 
u 2 (22/0961) 
E 2x0 10920 7 
[a] = E 0 0 -ጽ (2.1. 8) 
Оо а 
| 0B Wo (0 f 
à 91 0 
2 9. x = |6852 5: የ 
2 21. 0 
| Е | % B 





2.2 Transformation to Canonical Form 


Assuming in our linear transformation matrix [A] there is some 


eigenvector v associated with eigenvalue A. 


ji Mg 
0 0 NT di 
0 0 0 Ë -Л 


There comes the characteristic equation: 
ER ае 
The eigenvalues are: 


- 4.15 


2 
KS 


+ 4, 15 


x 


ru =-2 / 


Only А, апа ^ аге non-negative or unstable. 


where V & 0 


With these eigenvalues all distinct, one can always find a new set 


of state variables: 





related to x by the transformation 


such that the system of equations (2.1. 6) transforms to: 


ሜሽ О 
A 


ሃ = A cU (2.22 


with the 4 x 4 matrix [G] and vector D given later. 


From (2.2.1) 
x [e] ! y (2. 2. 3) 
Let 
511 812 813 814 
821 892 823 594 
[e] 71 > 
1] 


-1 
The transformation matrix [ገ] is not unique, but a convenient 
form for this problem can be found by four column vectors, all are 
eigenvectors associated with one eigenvalue respectively. 
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Define 


Vo УЗ 74) (2. 2. 4) 
- 0 MELA SS 4 (2. 2.5) 
But = NEU ሦር ; Na = + e 1. 8) appears as: 
1 0 0 
= 
0 07-24 
4 (2.1. 8) 
0 0 1 
0 0 ኣፈ 
= 1 
0 0 , 8 
_ 221 
e CE 521 
=0 
-A 
L | 531 


一 入 
人 


There are many possible solutions, one set of which can be: 





= 0 
541 
a 
ЕЁ 
21 Ar”, 
= -- Жа 
511 ^^ 
for i = 2, 3, 4 respectively (note 入 3 = 0), we can adopt: 
Bag © | 
= 0 
$32 
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5. 















2220 SAT 
843 ^ 0 
UN 
553 4 
ከጠ 
513 ` 9 
- ፪/ 
- а 
544 4 
mE 
> IRA 
524. (ACA4 Ax AQ) 
g = -入 1 
14 (Agra (A2-Aq) 
Then 
-入 2 ^j 0 те А) 
A2۸1 AGA Are) (^224) 
-MAz Ai^z 0 -Ла А 
A^ АЖАР CA,- Na )(A2-A4) 
a 9 9. 
еп : ^+ Ax 
0 0 0 ERE 
ይ 
Inverse of [G] is [c]: 
1 0 SIRE 
= Ay 4 21-74 
1 0 mr ^2 А4 
к С d ^ra (2.2.6) 
[°] | Ла и 
0 0 a з 
0 0 0 Ла 


ІП 





By (2. 1. 6) 
күүс, (2.1. 6) 


By transformation 


y 16] 
[Ie] “i = [al y + eu 
y - [ጣ[ዘጋ[ጣ ”፲ ፥ [ጣረ። (2. 2. 7) 


li () 
እን 
А; 


[6] (4] [a] > DES. 9 (2. 2. 8) 
Ла 
а D 
1 1-2 
d ange 
2 1 - , 
р = 
K Ра “| ) (2. 2. 9) 
d, ነ Да 
Finally, we get the Jordan Canonical form of the system: 
y СЕ 2: (2 2 10) 
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2.3 The Reduced System of Equations 

In the last section, the whole system is represented by equation 
(2.2.10). Now we have to show that for purpose of establishing a 
controller which assures stability of the equilibrium point at the origin, 
it is sufficient to consider the following reduced system. 

ነሃ ው diu i = 2,3 (29589 
without regard for the coordinates yi and Ул» which associated 
with negative eigenvalues. 

If this is true, i.e., a controller u = f (Yo, Уз) can be found 
which makes the system (2. 3.1) asymptotically stable for some region 
about the origin of the two dimensional space. By definition: 

У 7.9 as (%--> ОО 21 702 25 


Erom (2.3. 1) 


(у; - du) >0 as t— oo |. 
Then, certainly, as t—% 
t+ot 
ы (X-du)it 0  &*0 ,.,, 
ው 
But {At t4at 
Ғ-ҒА y, (tt50)- Y, (t) 4 1 dt 
Lim L| (y-diu)dt =Ltim | —— —- iat а 
т OF i TE ii £ 
t 
'Then ttal 
е 


|) 





This impliesu (t) —0 as t —-eoin the sense of its average 
value over any very small non-zero time interval. 

Now, we come back to those equations for 71 апа Уд: By 
means of Jordan canonical form, the state variables are express- 


ed by partial fraction as the following block diagram. 





4 
E ሃ 
ENS ! 
1 
а, 
ሃ 
SER 2 
u 
С 
d, 7 
TEN 4 


FIG. (2-3-1). State Variables in Jordan Canonical Form 
Let h; (t) be the impulse response for yj then 
Lim y. (9 = Lim |n (EISEN ml . = 1,4 
t>% J ር ጋ. 9 1 1 
Since 5 (t) approaches zero exponentially as t > oo 

(because of negative eigenvalue), and since as mentioned above, 
u (t) can be considered to approach zero under any integral sign, 
it follows that y; 一 > 0 аз 1-> oo, 
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ሽ 
D 
i 
" 








Thus the initial displacement of the states Yi and 7А һауе 
no effect on the region of stability of the complete system. There- 
fore, from now on, we only consider the reduced system described 


by equation (22391). 
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3. Control with an Analog Computer 
3.1 Selection ofa Controller 

A controller is needed to provide stability in the region of 
controllability, which means the largest region in the state space 
within which the system can be brought to the point of equilibrium 
at the origin with the constraint uU. The controller will be a 
function only of yo and Уз › but through the transformation y= [G]x 
it will generally involve all the state variables of the original 
system. Pontryagin's maximum principle determines an optimum 
u (t) which minimizes: 

T 

y, T) -J f (yo, y.) di 
with the constraint 

Jul <u 
and final states 

7, (ቸኃ ። 0 1 = 3, 3 

> is some positive cost function, different kinds of which 
lead to different kinds of controller. For the case of minimum- 
time controller m = 1. We define a new state variable: 

zu 
rg > | at 


It follows 


By adding this new state variable, our system gets: 
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(2፡41፡ 688) 


<< 
CD 
]| 
© 
። 
ፎ 


Pontryagin further defines the Hamiltonian function, maxi- 
mizing this H function with respect to u has the same effect as 
minimizing y Ct). 

In our case: 
4-2. P, y, =P, (Mht dzu)tPzurk, 

For maximization with respect to u, probably the bang- 
bang type controller (u = x U) is the best consideration. 
an 

u = U sgn [f] 


Usually f' can be derived by solving the following canonical 


equations. 
vi OP; 
б од 
P 22 


In this problem, it is hard to express them explicitly. In 
view of piecewise linear switching, the following guess seems 


reasonable, 
f = ከሃ - 
У #5 


where b is some positive constant. 


612 
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Therefore 

u = Usgn [by。 - yy] GENES 

If there is no control on cart position, that means zs no longer 
a state variable. Then by equations (2. 2.1) and (2. 2. 6) the 
uncoupled state Уз has no meaning. In this case the whole system 
(2. 3. 1) reduces to: 

уо = Aaya t do U (3.1.4) 

u = U sgn [ -yo) (SO) 
3,2 Realization of the Control 

Now, we get a minimum-time controller, which is (through 
the functions of yy and Уз) in terms of the original state variables, 
namely, 0i, Ө», E ¡ and Yo. Those state variables сап be gener- 
ated by two potentiometers and two techometers. We selecta 
Donner Analog Computer to sum them up and use a D. С. relay to 


generate sgn function. The real structure is shown as Fig. (3-2-1). 


3.2.1 No Control of Cart Position 





By (3.1.5) 
u = U sgn (- yo) (3 ТЕО) 
By ze 2. 6) 
я 1 1 Же 
о AA Ы 9 A24 52 


52፡0 0, - 0.138 5 (322 118) 


2. 
18 
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TT 
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"rolled 








By measuring those potentiometer and tachometers, the 
following data are obtained. 


08. ` 1 volt corresponds to 0. 087 radian 


5 1 volt corresponds to 40.4 radian/sec 


E 1 volt corresponds to 0. 36 meter/ sec 


Multiplying these factors, we get: 


u = U sgn [ 0. 087 Ө; + 9. 75 9, + 0. 05 É , J СІП 
where 
V fy 
Ú = —— = 0.24 ` (322441803) 
4 
The circuitry is built for the controller as shown in Fig. (3-2- 14), 
O5. 
0975 
1M 5 (v] 
Өс 
я 01M 
97 , 0.1/М 
10 
10M ARelay Coil 


5; 


0s connector Š 


FIG, (3-2-1-1). Controller Circuitry 
3.2.2 With Control of Cart Position 
[C SM) 
u = U sgn [bys у, | (3,186) 
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ETO. 2.6) 
ET 


T 


1 

= + — = N 
ye 2, 4 g,- 0.204 + 0.1028, (3.2.21) 
For s 

1 volt corresponds to 0. 05 meter 
Then 

= 0.010 L 08056 
y = A 


The following network is added to the connector in Fig. (3-2-1-1) 


52" сад 





To Connector 
of Fig (3-2-1-1) 


FIG. (3-2-2-1). Additional Network for Position Control 

If the value of b is not big enough, the cart has a tendency 
to settle itself a little bit off-center. For instance, taking b=0. 08 
(as shown in Fig. 3-2-2-1), the cart tends to settle itself about 0. 2 
meter from its starting position. If the polarity of the potentio- 
meter for 2 is reversed, the cart will tend to settle in the opposite 
side from the starting position. Anyhow, if b is big enough, the 
cart will come back to its starting position. This is shown in 


appendix Il-5, for b = 0.7. 
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4. Simulation by Graphs 
4,1 The Region of Controllability 

When we are playing with the physical cart, we can start the 
cart motion by just pushing the pendulum off its vertical position. 
Now, for the case of simulation with digital computer, we encounter 
the problem of deciding the initial condition of 91- In other words, 
we want to know the region of controllability which means the 
largest region in the state space from which the system still can 
be brought back to its equilibrium point. 

With a bang-bang type controller (u = * U), the trajectories 
consist of two segments corresponding to u = U andu = -U 
respectively. The origin of the two-dimensional space can always 
be reached by this switching of u. The region of controllability, 
then, can be defined as the set of points reachable by the trajec- 
tory starting at the origin (yy = 0, y3 = 0), and proceeding in 
reverse time(0St <- ©) with u alternately taking values of +U 
and -U. In our problem: 

Уә = Лә уд + йди 

(2095 1) 

Ya = dau 

` Those first order differential equations can be easily solved 
with solutions as: 


22 


159 
jg REED 


= exp[A,(t- to] (4.1.1) 
ከ. ከካ. 
E (4.1.2) 
ази ዕ 


From (4.1.2), Уз is undefinted as t — -oo, no matter what the 
initial value Уа0 is. Thus the region of controllability is unbounded 
in the y4 coordinate. Equation (4. 1. 1) shows directly that 
dou 

175 


d 
I*z| < 


— 











Уә as t— -0O, hence y, must be bounded by: 





Numerically we have 
|y5| -< 0. 45 (4. 119) 


By (2. 2. 6) 


2 1 2,124 
Y9 7 798 "X38 *3 ۸2-۸4 


If only 9 has non-zero initial value, it must be bound as: 


= 0.45 





|e (0 
This is the reason of using 0. 1 (radian) as the initial value of 
angle displacement. 
4.2 Analysis by Graphs 
4.2.1 No control on Cart Position, Equation (2.1.5) can be 


written as: 
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ነና е 
= 
© 
= 
© 
ርጋ 
ዞና 
pd 
O 


Xo 72 O 0 3.5 хо I7 

| — Х + ጃ u na L: 5) 
х. 0 0 0 1 ха 0 

x, 0 0 0 -г2 x, 9,8 


!ሸ ህብ...) СІРГЕ 

9 = 0.24 sgn ( x, + 0. ОЕ J (4.2. 1. 1) 

Now, this system can be simulated. For better ассигасу, ме 
use sin x, instead of xj. Graphs are shown in Appendix I. Same 
initial angle for all graphs. 

Fig. I-1 shows pendulum angle vs. time, with initial angle 
displacement 0.1 radian. It comes back very quickly. 

Fig. I-2 shows the pendulum changing rate vs. time. 

Fig. I-3 shows the cart position vs. time. 

Fig. I-4 shows the phase plane (8 vs. 0). The tracjectory 
does come back to the origin. It means stability. 

Fig. 1-5 shows ou ይ. Astime goes on, the position 
rate decreases linearly to zero. 

Fig. I-6 shows the control force vs. time. Itis the bang- 
bang type force; only the direction of the force is switched by the 
relay. 

4,2,2 With Control on Cart Position. 

In this case, the system equations are same as before, but 

the control force changes as: 


24 





DU 24 sgn| x,* 0. 241x9+(0. 102b+0. 138)x4+0. 204bxə J 


We simulate this system with the initial angle displacement 
as before (9 (o) = 0. 1 radian). Graphs are listed in Appendix II. 

Firstly with b = 0.1 for position control. 

Fig.II-1 shows pendulum angle vs. time. It is stable. 

Fig.II-2 shows the phase plane. The tracjectory goes to the 
origin. 

Fig. II-3 shows the cart position vs. time. It does not come 
back very quickly, because the amount of b is too small. 

Then take b = 0. 7 for position control. 

Fig. II-4 shows the stick angle vs. time. 

Fig. II-5 shows the position vs. time. The cart comes back 


to where it started very quickly. 
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5. Optimal Discrete-time Control 


9.1 The System in Discrete-time Form 


Recall the original system as follows: 


x= (A]x+cu (2. 1. 6) 
Consider, firstly, the equation without control force 
х = (А) х (5.150 
and assume a Taylor series 
(3.122) 


፳ (1) = ልፊ + ጴ11 + Aot® +...+A,t™ +... 


to be the solution of the above homogeneous differential equation. 


Then, set t = 0, one obtains: 


x (0) = A. 
Next, if (5. 1. 2) differentiated and then t set to zero, one obtains 
Ее. 
But, from (5.1.1) 
х (0) = [А] х(0) 


Then 


ew (А) х (0) 
If (5. 1. 2) is differentiated twice, and t set to zero, one 


obtaines: 
x (0) - 2A, 

Or 
2А, = Х (0) = (А) х (0) = (А) х (0) 


26 








So that 
m! 2 
Аз = lj[A] * x (0) 
Continuing this process, all the terms fig can be evaluated, 


and one obtains: 
сл па па Е 
xe «JG (а) GEE Маса Оаза ec) i o) (59123) 
n 2! m! Е 
By comparing it with the scalar expansion of eat, it is obvious 
to have a more compact form, like: 
x (t) = eAt x (0) (б ЗАП) 
In «my case, eAt isa 4x4 matrix. It is usually called funda- 
mental matrix and designated by: 


$ (t) Pe At 
Apparently 
J (=t) = A 077 (t ^ (8.1.4) 
Also 
x (t) - 0 (t) x (o) 
Then 
x (t) = (0 х() -(А)х(0 = (А) (9 x (o) 
So that 
| (t) = [AJ f (t) (5.1.5) 
In order to solve the equation (2. 1. 6), we try to find a particu- 


lar integral in the form of 
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x (t) = f (t) y (t) 
p 

By putting into (2. 1. 6), one obtains 

Ú t) y (O +f G) yQ =[A] J G) y GO + e u (t) 
ES I5) 

Ø (t) y (t) = e u(t) 

f -l 
y(t)=f # ወዩ #ጦ።#፻ 


Then 
£ ዓያ 
x, (0 = f(D % (2s (t)dT 
By (5. 1. 4) 
8 . 
x, 6 01 F-Yewr)dr (5. 1. 6) 


In evaluation of $ (t), we know the argument t represents the 
time interval between two instants. More conveniently, we can 
describe by: 

s LA MON) 

x (t3) - f (tg- t x (tg) < | (t3- to x (t) 

But X(6)= $(6-6G)$ (G-6) X (0g) 
рез) =$ (t3- to) Ba tp 
By this reason, (5. 1. 6) can be put into the more familiar form of 
a convolution integral. 
o (5.1.6.1) 


Then the general solution of (2. 1. 6) will be: 
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‚ 
mm 
м, Е 
mE БЕРІ 
p 





х (0 в (хо) + |  а-Әса(с)ат (о ШЕТ) 
о 
In case of discrete time, it turns out to be: 
x (k+1) = (DT) x (k) ዛማ (TW See me) we ንቀን) 
о 
Where 
D | 
DT = sampling time 
By noting of a constant control force through the interval DT, 


one obtains: 


x (k +1) = P (DT) x (1) + зо: | PT фот )сат (5.1.9) 
Ог 

x (k + 1) = # (0ፐ)- ፳ (k) +Á u(k) (5. 1. 10) 
Where 


DT 
42] # (Т -7) сат 
г 2 2 m m 
f ወ፲) ።[ ፲] +[ል] - 5፲ + [AJ (IDT) + (DT +... 
2! m / 
The computation of 2 (DT) and # (DT) would be very tedious, 
but, with the high speed computer, they can be obtained within 


seconds. The program to achieve this is shown in Appendix IV-3. 


(DT = 0.1 sec. ) 


Where 
-0.0817 
-1. 6067 
А = 
0.0458 
0. 8882 
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1. 0872 0.1029 0. 0000 0, 0166 


1. 7697 1. 0872 0. 0000 0.3268 
0.0000 0.0000 1. 0000 0. 0906 
0. 0000 0. 0000 0. 0000 0. 8187 


9.2 Evaluation of Control Force. 
In optimal discrete-time control, if we want to minimize 
the following cost function. 
J(n) z[: (к) Ох (к) е га? (к-1)) (5.2.1) 
The second term represents the amount of control force 
which can be allowed, arbitrarily r setto 1. The first term 
gives the choice of state variables which will be minimized. In 
my case, x, and X4 are those variables. So, Q becomes: 
1 Ж О 
elo“. 
By (S110) (5221) 


J( n) = fx(n-1+au(n- 1) Q |fx(n-1) +А u(n-1) Sen (eet + J(n-1) 


5871.) 
Qu (n-1) 


J(n-1) is independent of u(n-1), gives 


Minimizing with respect to u(n- 1), - 0 and noticing 


pa (n-1) f x- DA! дана Их (n-1) tau(n- 1| +2ru(n-1) = 0 
t 
Du - اود‎ _ x (n 1) (5. 2.7) 
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Define 


L 
at P-A QP (52245) 
1 A QAT 
Then 
t 
uln-1) = а х man) О, 
Substitute u(n-1) to (5. 2.1), one obtains: (5.2. 4) 


J (n) -fg (Их (а-2)+ди(а-2) + да “ y a 4| | + да i} 
+rf 于 et Oia Jt Q[ J+ ru(n-2) + I(n-2) 
Where ау = (а)! 
( JP[ f x (n-2) tou (а-2) | 
Further defining 
44 Оф +да ! (5, 2. 5) 
The first and third terms of (5. 2.4) can be combined as: 
[p ፲"ቻ+፤( 7". ዛ|#.[ ] ፅጪቫ }}+ ]'9፪ ] 
Сои + одана, а ОД за дола JH JO 
A JUPO aa O A O AN 
MAS O a) AO 
o 
=( Jw Q# + QC J 
Now, (5. 2. 4) becomes: 
та) Г | оў +оо ЕИО СЕ 
[1*9 Qu, +0+rajaf)[ )+ ruó(n-2) + J(n-2) 
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Define 


D Ж t 
Br Mo S ааа 220 


m) 
2 u(n-2) 


- 0 gives: 


[ 1 ра *A 1 J + 2 ru(n-2) = 0 


(ES (n-2) gt TU (п-2)4')р д +atp L 1 + 2 ru (n-2) = 0 


[4 
u (n-2) = - ай. x (n-2) СО ОТ) 
АР АТТ Бо 
Define 
t p, f 

a, = Шр (5-22 
Then 

9(8-2) = 83. ፳ (8-2) (5 2100) 
Define Po D С) 


(5. 2. 3) becomes: 


t 
NM A p f (5.2 
1 Ai pA in 
O 
Continuing one more stage, and set gJtn)_ - 0, one obtains: 
229(-3) 
Е t 
y = f 4а; 
2 t t 
Pa Yo P1jo* Q* rag 
AD, Y 


t = MALA 
an А род г 


u(n-3) = a, x (n-3) 
Continuing on the same procedure, one can expect the 
following general forms. 
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ат t 
Pk “Yk Pk-1 Ret Qtr aya, 


t 
t AP f _ 


ыз. መበ 

u(n-k) = a х (п-К) 

We note that u(n-k) depends solely on those present states 
х (п-к). Thus makes clear Bellman's ’'Principle of Optimality", 
which states: ''An optimal policy has the property that whatever 
the initial state and the initial control input rector are, the 
remaining control input rectors must constitute an optimal policy 
with regard to the state resulting from the first control signal. " 

When the number of stages gets very large, a converges 
to some final set of values. Thus exists some fixed values for 
the feedback compensation for all values of time and state variables. 

With the help of CDC 1604 computer, 200 stages are accom- 
plished. The program is shown in Appendix IV-4. Finally, one 
obtains, for all sampling time: 

2, 4846 


0, 5990 
u(t) = x (t) (5. 2. 9) 


0. 0074 





0, 3448 | 
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9.9 Simulation by Integration 


With the control force shown in (5. 2. 9), we can simulate 
the system by calling a subcontine for integration. The control 
force is calculated during every sampling instant and hold con- 
stant until the next sampling instant. 

The FORTRAN program achieving this purpose is shown 
in Appendix IV-5. Initial angle displacement is same as before, 
x (1) = 0.1 radian. Sampling time DT = 0.1 second. 

The graphs are collected in Appendix III. 

Fig. III-1 shows stick angle vs. time. The angle comes 
back after about 40 stages. 

Fig. III-2 shows phase plane. 

Fig. III-3 shows control force vs. time. 

Fig. III-4 shows cart position vs. time. It gets back to the 


starting position comparatively slowly. 
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APPENDIX I 
Graphs for No Control of Cart Position 
stick Angle vs. Time Without Position Control 
Stick Angle Changing Rate vs. Time Without Position Control 
Postion vs. Time Without Position Control 
Phase Plane Without Position Control 
Position Changing Rate vs. Position Without Position Control 


Control Force vs. Time Without Position Control 
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A-SCALE = 1,88+28 UNITS INCH. NDS 
M-SCALK = 92.20€ -02 UNITS/INCH. 


WAN NO POSITION CONTROL 
RUN 1 STICK ANGLE US T 


X-Scale: 1.0 Second/Unit 
Y-Scale: 0.02 Radian/Unit 


FIG. (I-1). Stick Angle Vs, -Time Without Position Control 
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XSCALE = jQQE+22 UNITS~ INCH. 
M-SCAK - laM -8 UNITS INCA 


WAN NO POSITION CONTROL 
RUN 1 ANGLE RATE US T 


X Scale: 1.0 Second/Unit 
Y Scale: 0.1 Radian per Second/Unit 


FIG. -(I-2). Stick Angle Changing Rate Vs. Time Without Position Contre: 


vf 





935 


234 


225 


215 


918 


225 





X-SCALR = 1.00€£+208 UNJTS/ INCH 
M-SCALK = SARE-RR UN] T57 INCA 


WAN NO POSITION CONTROL 
RUN 1 POSITION US T 


X Scale: 1.0 Second/Unit 
Y Scale: 0.05 Meter/Unit 


FIG, (1-3). Position Vs. Time Without Position Control 
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XSCALE = 2.808€ -92 UNITS INCH 
'1-30ዲቪ = 1,80€-81 UNJTS/ INCH. 


WAN 
KUN 1 


X Scale: 
Y Scale: 


NO POSITION CONTROL 
PHASE PLANE 


0.02 Radian/Unit 
0.1 Radian per Second/Unit 


FIG. (-4). Phase Plane Without Position Control 
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A-SCALE = .5.ደደ -ደ2 UNITS INCH. 
M-SCALE = ፉጌ.ይጄ -ል1 ዚዮህ? 2።፲ኮርዮኒ 


WAN NO POSITION CONTROL 
RUN 2 ZDOT US 7 


X Scale: 0.05 Meter/Unit 
Y Scale: 0. 1 Meter per Second/Unit 


FIG. (1-5). Position Chdngihg Raté Vs. Position Without Position Control 
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X-SCALE = ፉ1ይደደ4ደፅ UNITS INCH 
ተ-53ርጧ.ጀ = JAQE-21 UNITS“ INGA 


WAN NO POSITION CONTROL 
RUN 9 PASA 


Note : Control force has no dimension. 


X Scale: 1.0 Second/Unit 
Y Scale: 0. 1 Umt Unit 


FIG, (1-6) Control Force Vs. Time Without Position Control 
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APPENDIX II 
Graphs for Control of Cart Position 
Stick Angle Vs. Time With Less Position Control (b = 0.1) 
Phase Plane With Less Position Control (b = 0.1) 
Position Vs. Time With Less Position Control (b = 0.1) 
Stick Angle Vs. Time With More Position Control (b = 0.7) 


Position Vs. Time With More Position Control (b = 0.7) 
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ASCAR = 3,2848 UNJTS/ INCH. 
“SCALE = 92.80€ -92 UNJTI7 INCH, 


WAN WITH POSITION CONTROL 
RUN 1 STICK ANGLE US T 


X Seale: 1.0 Second/ Unit 
Y Scale: 0.02 Radian/ Unit 


FIG. (II-1). Stick Angle Vs. Time With Less Position Control (b=0. 1) 
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ХАЛ * '..ዳደደ -ይኃ ዚዘህ ፲ ኃ።፲ከርቡኒ 
M-SCALR = J,QQ-21 UNITS INCH. 


WAN WITH POSITION CONTROL 
RUN 1 PHASE PLANE 


X Scale: 0.02 Radian/Unit 
Y Scale: 0.1 Radian per Second/Unit 


FIG. (II-2). Phase Plane With Less Position Control (b=0.1) 
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A-SCARK = LARE+BB UNITS/INCH 
M-SCALIE = S,2QE-22 UNITS/INCH 


WAN WITH POSITION CONTROL 
RUN 1 POSITION US T 


X Scale: 1.0 Second/Unit 
Y Scale: 0.05 Meter/Unit 


FIG. (I-3). Position Vs. Time With Less Position Control (b = 0.1) 
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K-SCAR = 1.80€488 UNTS/ INCH, 
М-5СА К = 2.09 -42 UNIT S7 INCH. 


WAN WITH POSITION CONTROL 
RUN 2 STICK ANGLE US T 


X Seale: 1.0 Second/ Unit 
Y Scale: 0.02 Radian/Unit 


FIG. (II-4). Stick Angle Vs. Time With More Position Control (bz0.']) 
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A-SCAR = 1.RABE+RAB UNITS/INCH 
M-SCALK = 5,ARE-22 UNITS INCH. 


WAN WITH POSITION CONTROL 
KUN 2 POSITION US T 


X Scale: 1.0 Second/Unit 
Y Scale: 0.05 Meter/Unit 


FIG. (I-5). Position Vs. Time With More Position Control (b - 0. 7) 
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APPENDIX III 
Graphs for Discrete-time Control 
Stick Angle Vs. Time for Discrete-time Control 
Phase Plane for Discrete-time Control 
Control Force Vs. Time for Discrete-time Control 


Position Vs. Time for Discrete-time Control 
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X-SCALE = 1,00£+20 UNITI/ INCH, 
MM SCALE = 2.80 -82 UNITS INCH 


WAN DISCRETE TIME CONTROL 
RUN 1 STICK ANGLE US T 


X Scale: 1.0 Second/Unit 
Y Seale: 0.02 Radian/ Unit 


FIG, (III-1). Stick Angle Vs. Time for Discrete-time Control 
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X-SCALK = 2,.QQ6-22 UNITS/INCH 
1-SCAR = 9.04-22 UNITS/INCH 


WAN DISCRETE TIME CONTROL 
RUN 1 РАЗВЕ 


X Scale: 0.02 Radian/Unit 
Y Schle: 0.05 Radian per Second/Unit 


FIG. (III-2), Phase Plane for Discrete-time Control 
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X-SCAR = 10-00 UNITSZINCH. 
M4 SCALE + SALE -22 UNJTS/ INCH, 


WAN DISCRETE TIME CONTROL 
RUN 1 | 


Note: Control force has no dimension. 


X Seale: 1.0 Second/ Unit 
Y Scale: 0.05 Unit/Unit 


FIG. (III-3). Control Force Vs. Time for Discrete-time Control 
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A-SCALE 7 2,00£ 481 UNITS/INCH 
М-5с4 К = 5.АЖ-02 ІМІТ5/ІМСН. 


WAN DISCRETE TIME CONTROL 
RUN 2 POSITION US T 


X Scale: 20.0 Second/Unit 
Y Scale: 0.05 Meter/Unit 


БІС. (1-4). Position Vs. Time for Discrete-time Control 


52 





APPENDIX IV 
FORTRAN Programs 
FORTRAN Program for Simulation without Position Control 
FORTRAN Program for Simulation with Position Control 
FORTRAN Program for Calculating ф and A Matrix 
FORTRAN Program for Calculating Discrete-time Control Force 


FORTRAN Program for Simulation of Discrete-time Control 
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С 


PROGRAM ВКООМ ' 

NO CART POSITION CONTROL 

DIMENSION X(30)»XDOT(30)»>»C (15) ` 
Cc(10)=1+0 

CALL INTEG1 (T»X$»XDOT»C) 

XDOT(1)2zX(2) 
XDOT(2) 17. 2*SINF(XC1) ) *3e95*X(4)-176e2*U 
XDOT(3)=X(4) 

XDOT(4)=-=2.0*X(4)+9.8*U 

U=z0.+24*SIGNF(lo »X(1)+C(1)%*Xx(2)+C(2)*X(64)) 
X(5)=U 

со то 1 

ЕКО 

END 
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PROGRAM BROOM 

POSITION CONTROL 

RUN NO. TWO FOR MORE Y3 FEED BACK 
DIMENSION X(30)»XDOT(30)»C(15) 

С(19)1. 

CALL INTEG1 (Т»Х»ХООТ»С) | 
XDOT(1)=X(2) 

XDOT(2)217.2*SINF(X(1))434 5*X(4)—1742*U 
XDOT(3) =X (4) 

XDOT(4)=-2.#X(4)+9.8*U 

0-0. 24#5106ዘ፻(1. )፤(1)+ር(1)*፳(2)+(ር(3)*#0. 102121) 
1%Х(4)%С(3)%0.204%Х(2)) 

ЕО 

GO TO 1 

END 

END 





e. JOBO250F + МАМ 


С 
E 


1003 
P 


аз 


400 


500 


401 


PRUGRAMSPHIDEL 

COMPUTING PHI AND DELTA MATRIX 

DT=SAMPLING TIME | 

DIMENSION A( 12912) »PHI (12912) »TERM(12912) »WORM( 12512) 


AE ИСА АТУ ТЕМ САА ር...) 


N24 

рТ-0.1 

ТМ= 0.0 

READ1>((A(IR>IC)>IC=19N)>IR=19N) 

FORMAT ((4F20.8)) 

ҺЕАр 1>(С(1)»91-1»М) 

PRINT 399sDTs((A(IR>IC)s»IC=1»N)»IR=]1;>N) 
FORMAT(///1X»3HDT=»1F5.3///»1X» 7HA(I»J)=/,tt4F10.2))) 
PRINT 3991 (C(1)s1=]p N) 
ЕОКМАТ(///1Х»5НС(І)-/(4Ғ10.2)) 


DO 400 IR=15;N 

DO 400 ІС-1»М 

TERM(IR»IC)=0.0 - 
WORM(IR>IC)=0.0 

TERM( IRs IR) =1¢0 

TELMC IRs IC) =TERM(C IRs 1C) *DT 

DELP(IR»IC)=TELM(IR> IC) 


ОЕМ ТК,ТС) 20.60 


РЕГ (ТК) 50, 
РН! (ТК „ТС ) -ТЕКМ ( ТК, 1 С) 


TM=1.0+TM 

DO 500 1=1 5 ከ 

DO 500 IC=1 sN 

DO 500 ӘМт-1»М 
DELM(LIRSIC)ZDELM(IR$ IC) -TELMCIRS JN) XACAJNS ICI XDT/CTMT1«0) 
WORM(IR$IC)ZTERM(IRS JN) *ACUJNS IC) XOT/TM-*-WORM(CIR9IC) 
DO 401 ІКг1»М 

DO 401 IC=15N 

TERM(IRsIC) =WORM( IRs IC) 
TELMCIRs IC) =DELM( IRs IC) 
DELP(—(IRsIC)=DELP(”(IR,+IC)+TELM(IR + IC) 
РНТ (ТК. ,ТС) +РНТ ( ТК,ТС ) + ТЕКМ ( ТК, ТС) 
ОЕСМ (ТК ,ТС) OS 

WORM(IR»IC)=0.0 

АВС-0.0 

DO 2 I=1l+N 

DO 2 J=1 +N 

АА<ТЕКМ (Т, 7) 

AB=ABSF (AA) 

ТЕ (АВС-АВ) 35392 

FIND BIGGEST VALUE 

ABC=AB 

CONTINUE 

ГЕ СО. 000000005 =АВ ЈЕ 5. 

GO ТО 4 

PRINT 502»((PHI(IR»IC)»IC=1»N) »IR=15>»N) 





2 са ОНОК И АТ ХЕЛЕНЕ Oy 


рО 600 ТЕ, М 
DO 600 К=1, 
DO 600 J=19N 

600 DEL{(1)=DEL(1)+PHI( 19d) ¥DELP(JsK)*C(K) 
PRINT 503 (DEL(I)sIz15N) 

503 FORMAT(///1X>THDEL(I)=//(4F15.9)) 


END 
END 
DT=0.10 
A(I，J) = ү 
.00 1.00 90.0 e O00 
17.20 00 e 00 2.50 
4900 . 00 .00 1.00 
.00 00 .00 -2.00 
ር(1)= | 
„ОО -17.20 00 9.80 
PHI( I,J) = 
T08 123206 03102891421 0. 000000000 
le 769732445 Ше 56 Ое. 000000000 
0000000000 0.000000000 1.00000000 
0.000000000 0.000000000 000000000 
DEL(I)= | 
-.081750092 -1.606739468 0.045890345 
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(15(]115533.1700 
О.З25 855 
«090634623 
«819720959 
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ee JOBO2Z50F SWAN 


С 
1 
2 
3 
44 
45 
46 
47 
5 
С 
С 
6 
" 
8 
10 
C 
C 
13 
C 
С 


PROGRAM ОРТСОМ 
MINIMUM Y=SUMIXT(K )*Q%X(K)+4RXUX%2) 
DIMENSION PHI (494) sPSI(494) »9P (494) sDEL(4) sAT (2094) > 


16GM(494) 9Q(494) 9FM(4) sEM(4) 


READ l1»N»NM,NPRINT 

FORMAT (8110) 

КЕАО2,К „ОТ 

READ2» ((Q(lI>»J)»J=1>N)>I=19N) 
READ2s( (PHI (IsJ),J=l9N) sI=19N) 
READ2>(DEL(1)>I=19N) 
FORMAT ( (4Е 16,9) ) 

PRINT 39NsNMsNPRINT 
FORMAT(///(8110)) 

PRINT &&s R DT 

FORMAT (//1X92HR=>F5.e2»3X>3HDT=»F5e2) 
PRINT 45>((Q(1l>»J)»J=19N)>1=19N) 
FORMAT (//1Xs7HQ(UI 9J) =7 ((4F1669) )) 
PRINT 465 ((PHI(CIsJ) sJ=l9N)sT=19N) 
FORMAT (//1Xs9HPHI(IsJ)=/((4F1629) ) ) 2 
PRINT 475 (DEL(I)sT=15N) 

FORMAT (//1X%> THDEL(C 1 ) =7 (4F lowed)? 

ОО бІс1»М 

DOS .ህ=1 5 እዚ 

GM(1»J)=0.0 

ЕМ(І)-0.0 

FM(1)=0.0 

P(1»J)=Q(I»J) 

Pol Clie) )=0 60 


CRECUEATE TATEN 

DO 22 KK=1>NPRINT 

ОО 20 Қ-1»ММ 

ОЕМ= 0..0 

006 I=1>N 

DO 6 J=19N 
EM(1)=EM(1)+DEL(J)*P(J9I1) 
DO 8 I=15>N 

DO 7 J=lsN 
FM(I)=FM(I)+EM(J)x*P Nn C Jo) 
DEN=DEN+EM(I)*DEL (I) 
DEN=ZDENZR 

DO 10 I=1sN 
AT(KsI)=FM(I)/DEN 
FM(1)=0.0 

EM(1)=0.0 


CALCULATE PSI(Ks+ 159900 

ОО 13 Іг1»М 

DO 13 JslsN 
PSI(IsJ)2PHI(IS;J)TDELCI)XATCR9 2) 


CALCULATE Р(К „То. 7) 
DO 16 [=l>5N 





ВО 168 1-51. 
DO 15 L=l+N 
DO 15 M=19N 

15 GM(I>9J)=GM(I>J)+PSI(L>I)*P(LIM)PSI(M>J)) 
Р(1»2/)-ОМ(1»>.)) «ҰВХАТ(К91)%АТ(Қ».) 

С FOR TERMINAL CONTROL OMIT Q(I15J) 
16 GM(1I5J)=0.0 
20 CONTINUE 


C р 
C PRINT 
PRINT 12»KK» (AT(INM»J)»J=15»N) 
12 FORMAT(//1X»3HAT(»12»2H)=/(4F16.9)) 
22 CONTINUE 
END 
END 
АТ(10)- 


2.484604417 «2990829066 «007459609 «344834834 
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ee JOBO250F »WAN 
PROGRAM BROOM 


C DINE AL DI SEGRETE ПО ስር ን ERDE 

ር РОБ СОТ МРЕТСТТУ, UÜSEFAT INT IN ASI ST Gr PE › ው 51 
DEMENSTON X130)5XDOT 120 1357) 
етет 


ПЕ САО ТЕСИ 4 1sXsXDOT sC) 
DOTA EZ) 
АВӘОТТОТЕТТ»2ХӘТИР(ХГІ)) 35 ZITO 
XDOT(3)=X(4) 
XDOT (4) =-2.*X(4)+9.8%*U 


ОТ-0.1 
E DT = HOLD TINE 
КЕКСИ 27 
С ABOVE STATEMENT AS A ZERO ORDER HOLD . 
2 U=C(1) ¥X(1)4C02)¥X(2)4+C( 3) 4X13 )40(4) or 
X( 5) =UÛ 
iS SD 
SO OD 
END 
END 
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